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ABSTRACT 

Galactic spiral shocks are dominant morphological features and believed to be responsible for substructure 
formation within spiral arms in disk galaxies. They can also contribute a substantial amount of kinetic energy 
to the interstellar gas by tapping the (differential) rotational motion. We use numerical hydrodynamic sim- 
ulations to investigate dynamics and structure of spiral shocks with thermal instability in vertically stratified 
galactic disks, focusing on environmental conditions (of heating and the galactic potential) similar to the Solar 
neighborhood. We initially consider an isothermal disk in vertical hydrostatic equilibrium and let it evolve 
subject to interstellar cooling and heating as well as a stellar spiral potential. Due to thermal instability, a 
disk with surface density Sq ^ 6.7 Mq pc~^ rapidly turns to a thin dense slab near the midplane sandwiched 
between layers of rarefied gas. The imposed spiral potential leads to a vertically curved shock that exhibits 
strong flapping motions in the plane perpendicular to the arm. The overall flow structure at saturation is com- 
prised of arm, postshock expansion zone, and interarm regions that occupy typically 10%, 20%, and 70% of 
the arm-to-arm distance, in which the gas resides for 15%, 30%, and 55% of the arm-to-arm crossing time, 
respectively. The flows are characterized by transitions from rarefied to dense phases at the shock and from 
dense to rarefied phases in the postshock expansion zone, although gas with too-large postshock-density does 
not undergo this return phase transition, instead forming dense condensations. If self-gravity is omitted, the 
shock flapping drives random motions in the gas, but only up to ^ 2-3 kms"' in the in-plane direction and less 
than 2 kms"' in the vertical direction. Time-averaged shock profiles show that the spiral arms in stratified disks 
are broader and less dense compared to those in unstratified models, and that the vertical density distribution is 
overall consistent with local effective hydrostatic equilibrium. Inclusion of self-gravity increases the dense gas 
fraction by a factor ^ 2 and raises the in-plane velocity dispersion to ~' 5-7 kms"'. When the disks are massive 
enough, with So > 5 M© pc"^, self-gravity promotes formation of bound clouds that repeatedly collide with 
each other in the arm and break up in the postshock expansion zone. 

Subject headings: galaxies: ISM — instabilities — ISM: kinematics and dynamics — methods: numerical — 
stars: formation 



1. INTRODUCTION 

Spiral arms in disk galaxies are regions of ongoing active 
star formation, sharply outlined by bright young star com- 
plexes. They usually span the entire optical disks and some- 
times extend even to outer gaseous disks (e.g.. Dickey et al. 
1990; Boomsma et al. 2008; Berlin & Amorisco 2009 and 
references therein). Such large-scale spiral patterns may be 
the manifestation of spiral density waves which propagate 
with a constant pattern speed through stellar disks (Lin & Shu 
1964, 1966), or may be transient features driven, for example, 
by tidal interactions with companion galaxies (e.g., Toomre 
& Toomre 1972; Hernquist 1990; Salo & Laurikainen 2000; 
Oh et al. 2008; Dobbs et al. 2010). Regardless of the ori- 
gin of spiral features, it is widely accepted that the interstel- 
lar medium (ISM) is strongly compressed when it encounters 
stellar arms, forming narrow dust lanes in optical images (e.g., 
Elmegreen et al. 2006; Shetty et al. 2007). The densest parts 
of the shocked layers subsequently undergo gravitational col- 
lapse and produce downstream secondary structures, such as 
OB star complexes and giant H II regions (e.g., Baade 1963; 
Elmegreen & Elmegreen 1983), filamentary gaseous spurs 
(also referred to as "feathers;" e.g., Scoville et al. 2001; Ken- 
nicutt 2004; Willner et al. 2004; La Vigne et al. 2006; Gordon 
2007; Corder et al. 2008), and giant molecular associations or 
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atomic superclouds (e.g., Elmegreen & Elmegreen 1983; Vo- 
gel, Kulkarni, & Scoville 1988; Rand & Kulkarni 1990; Koda 
et al. 2009). 

Since shock compression within the arms is the first step 
towards star formation in grand-design spiral galaxies, un- 
derstanding structural and dynamical evolution of these gas 
flows is essential to a host of fundamental problems, such 
as global star formation rates, the nature of the Hubble se- 
quence, galaxy evolution, etc. Since the pioneering work 
of Roberts (1969) who obtained one-dimensional stationary 
shock profiles, there have been numerous studies of the struc- 
ture of galactic spiral shocks under the simplifying assump- 
tion that the gas remains isothermal (e.g.. Woodward 1975, 
1976; Lubow et al. 1986; Martos & Cox 1998; Kim & Os- 
ti-iker 2002; Gomez & Cox 2002, 2004; Wada & Koda 2004; 
Boley & Durisen 2006; Kim & Ostriker 2006; Dobbs & Bon- 
nell 2006). In paiticular. Woodward (1975) and Kim & Os- 
triker (2002) showed that the one-dimensional shock profiles 
found by Roberts (1969) represent stable equilibria when the 
fluid quantities are allowed to vary only in the direction per- 
pendicular to the arms. The growth of axisymmetric self- 
gravitating modes is limited by postshock expansion (Lubow 
et al. 1986). 

When the direction parallel to the arm is included in mod- 
els, on the other hand, isothermal spiral shocks in two dimen- 
sions are prone to various kinds of instabilities. Balbus (1988) 
showed that when self-gravity is included, spiral arms are un- 
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stable to transient swing instability. When magnetic fields are 
included, spiral arms are subject to magneto-Jeans instability, 
in which embedded parallel magnetic fields that exchange an- 
gular momentum limit the stabilizing effect of galaxy rotation, 
encouraging non-axisymmetric perturbations to grow into gi- 
ant clouds and other arm substructures (Elmegreen 1994; Kim 
& Ostriker 2002, 2006; Shetty & Ostriker 2006). Wada & 
Koda (2004) showed that spiral shocks in two-dimensional 
thin disk models are unstable to a vorticity-generating wiggle 
instability and develop arm substructures (see also Dobbs & 
Bonnell 2006), although these in-plane modes appear to be 
suppressed by embedded magnetic fields (Shetty & Ostriker 
2006; Dobbs & Price 2008) or by vertical shear and mixing 
when all three dimensions are included in models (Kim & 
Ostriker 2006). 

While steady in-plane shock solutions are subject to insta- 
bility, shock models that include the vertical degree of free- 
dom do not even have steady solutions. Instead, the shock 
front in vertically stratified disks moves back and forth rel- 
ative to the mean position (Kim & Ostriker 2006). These 
shock "flapping" motions arise mainly because the vertical 
oscillation period of the gas is, in general, incommensurate 
with the arm-to-arm crossing time, so that the gas streamlines 
are not closed. Kim, Kim, & Ostriker (2006, hereafter Pa- 
per I) showed that the shock flapping is able to feed random 
gas motions on the scale of disk scale height that persist de- 
spite dissipation. The induced gas velocity dispersions reach 
a sonic level, suggesting that spiral shocks may be a con- 
siderable source of the ISM turbulence. Motions driven by 
shock flapping motions destroy any coherent vortical struc- 
tures that would otherwise grow near the spiral shocks, sup- 
pressing the wiggle instability. Since gravity is a long-range 
force and insensitive to small scale density structure, how- 
ever, magneto- Jeans instabilities still grow within the arms in 
three-dimensional disk models, in spite of non-steady motions 
induced by shock flapping (Kim & Ostriker 2006). 

Phase transitions caused by thermal instability (TI) create 
a multiphase ISM, with important consequences for galac- 
tic star formation. In the classical picture of the ISM, TI 
changes an otherwise uniform ISM into warm rarefied ma- 
terial and cold dense clouds in a rough pressure balance (e.g.. 
Field 1965; Field et al. 1969; Meerson 1996; Heiles 2001; 
Wolfire et al. 2003), while there also exists significant mass 
in the thermally-unstable temperature range (e.g., Heiles & 
Troland 2003). Supernova blast waves create an additional, 
hot component that is organized into bubbles or cavities (Cox 
& Smith 1974; McKee & Ostriker 1977), although the total 
mass contained in the hot phase is much smaller than that 
in the cold and warm phases (e.g.. Cox 2005). Cold atomic 
clouds transform to molecular clouds if their volume den- 
sity (to produce molecules fast enough) and column density 
(to self-shield molecules against photodissociation) are suffi- 
ciently high (e.g., Elmegreen 1993; Draine & Bertoldi 1996), 
as in, for instance, shocks in turbulent flows (Glover & Mac 
Low 2007; Glover et al. 2010). That the star-forming molec- 
ular clouds strongly correlate with spiral arms (e.g., Stark 
1979; Solomon et al. 1985; Kenney 1997; Zimmer et al. 2004; 
Shetty et al. 2007) suggests that spiral shocks too should trig- 
ger phase transitions from warm and diffuse to cold and dense 
conditions. 

Effects of TI on spiral shocks were first studied by Shu et al. 
(1973), who calculated one-dimensional shock profiles con- 
sisting of two stable phases in equilibrium. Although they 
aUowed for phase transitions, they assumed instantaneous 



thermal equilibrium, which precluded the existence of transi- 
tory thermally-unstable gas in their calculations. Using direct 
time-dependent numerical simulations including ISM heating 
and cooling, Tubbs (1980) and Marochnik et al. (1983) found 
that spiral shocks trigger phase transitions if the initial density 
is large enough. Because of strong numerical diffusion asso- 
ciated with insufficient resolution, however, they were unable 
to capture TI in the postshock transition zone, which is the 
thermally unstable regime. 

In Kim, Kim, & Ostriker (2008, hereafter Paper II), we used 
high-resolution one-dimensional simulations to study dynam- 
ical and thermodynamical evolution of gas flows across spiral 
arms with ISM heating, cooling, and thermal conduction. We 
found that even with TI, a quasi-steady state develops with 
the foUowing recurring cycle: both warm and cold phases in 
the interarm region are shocked and immediately transform to 
denser cold gas in the arm, which subsequently evolves to be 
Tl-unstable due to postshock expansion in a transition zone, 
and separates back into warm and cold phases. For a model 
with the initial number density of 2 cm"^, the gas stays in the 
arm, transition, and interarm zones for 14%, 22%, and 64% 
of the arm-to-arm crossing time, respectively. The gas mass 
in the Tl-unstable temperature range was ~ 25 - 30% of the 
total, and mostly located in the transition zone, suggesting 
that postshock expanding flows are important for producing 
intermediate-temperature gas. Paper II also found that TI in 
association with one-dimensional spiral shocks can drive ran- 
dom gas motions at ^ 1 .5 kms"' in the interarm and transition 
zones; these velocities are ^ 5-7 times larger than those from 
pure TI alone (e.g., Kritsuk & Norman 2002; Piontek & Os- 
triker 2004). 

In this paper, we extend the one-dimensional models of Pa- 
per II into two dimensions, in order to study the effect of 
vertical disk stratification on the dynamics and structure of 
multiphase galactic spiral shocks. The current work also ex- 
tends the vertically-stratified isothermal models considered in 
Paper I to include the effects of the ISM heating and cool- 
ing. Our key objective is to find how the flapping motions of 
spiral shocks inherent in stratified disks interact with multi- 
phase gas produced by TI, to change the shock structure and 
drive random gas motions in each phase. We also study the 
internal properties of clouds that form in self-gravitating mod- 
els. Although Dobbs & Bonnefl (2007, 2008) and Dobbs 
& Price (2008) ran SPH simulations to study shock struc- 
ture and cloud formation in three dimensions, they used 
pre-determined cold and warm particles and did not allow 
the transitions between them. Dobbs et al. (2008) included 
ISM heating and cooling in the energy equation and thus 
handled TI self-consistently, focusing on the formation of 
molecular clouds in spiral shocks. Using grid-based three- 
dimensional simulations, Wada (2008) studied dynamics of 
galactic gas flows under the influence of self-gravity, a spiral- 
arm potential, radiative cooling, star formation, and energy 
feedback from supernova explosions. Although these three- 
dimensional global models are of course more reaUstic, our 
local models are useful for studying the detailed dynamics of 
spiral shocks at high spatial resolution, and allow us to isolate 
each effect of the physical processes involved. 

This paper is organized as follows. In §2, we describe the 
basic equations we solve and specify the model parameters. 
In §3, we present the results of one-dimensional solutions for 
vertical disk equilibria including heating and cooling, also 
providing simple analytic expressions for the mass fractions 
and scale heights. In §4, we present the overaU evolution, 
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structure, and statistical properties of spiral shock flows with 
TI in stratified disks, based on the results of two-dimensional 
simulations. The effect of self-gravity is discussed in §5. In 
§6, we summarize our results and discuss their implications. 

2. NUMERICAL METHODS 

The local formulation used in the present study is similar to 
that in Papers I and II. In this section, we explain our numeri- 
cal method and model parameters. 

2.1. Basic Equations 

We consider a local region centered on a spiral arm that 

is assumed to be tightly wound with a pitch angle sin ; <Si 1 
and rotating at a constant pattern speed Vtp. We set up 
a local Cartesian frame, centered at the position {R.(p,z) = 
iRo,^pt,Q), that corotates with the spiral arm. The x- and 
y-axes of the local frame are aligned in the midplane par- 
allel and perpendicular to the local arm segment, while z- 
axis points the direction perpendicular to the galactic plane 
(Roberts 1969; Paper 1). Our simulation domain is a two- 
dimensional rectangular region with size L;^ x in the x-z 
plane (hereafter XZ plane). We assume that all physical vari- 
ables are independent of y (quasi-axisymmetric), while al- 
lowing nonzero velocity in the y-direction in order to han- 
dle epicycle motions associated with galactic rotation self- 
consistently. 

In this local frame, the galactic differential rotation is trans- 
lated into the background velocity as 

vo = Roi^o - ^p) sin /x+ [Roi^o -^p)- qo^oxjf, (1) 

where fio = ^(Ro) and = -idlnn/d\nR)\R^ denotes the lo- 
cal shear rate in the absence of the spiral potential (Kim & 
Ostriker 2002, 2006). Under the local approximation (i.e., 
Lx ^ Ro and |v| <C Ro^o), the equations of hydrodynamics 
expanded in the local frame are 



dp 
dt 



+ V-(pv7-) = 0, 



(2) 



d\T 1 

— - + Vr • VVr = VP - ^of^OVO;cy - 2f]o X V - V($, + $ext), 

at p 

(3) 

-— +Vr-Ve= '—PV-yr-pC, (4) 

at 7-1 



V^^, = 4nGp, 



(5) 



(e.g., Roberts 1969; Shu et al. 1973; Kim & Ostriker 2006), 
where vj^ = vq + v is the total velocity in the local frame, 
is the self-gravitational potential of the gas, $ext is the exter- 
nal stellar potential, and pC{p, T) is the net cooling function. 
Other symbols have their usual meanings. We adopt an ideal 
gas law P = (7- l)e with 7 = 5/3. 

The external stellar potential $ext varies in both the x- and 
z-directions and is separable into two parts as 



$ext = 2nGp^z + $sp cos — 



(6) 



where is the unperturbed midplane stellar density, $sp is 
the amplitude of the spiral potential, and = IttRq sin i/m is 
the arm-to-arm separation for an m-armed spiral. The first 
term in equation (6) adopts a constant density for the stel- 
lar disk vertically; this is a reasonable assumption given that 
most of the gas is located within one stellar scale height of 



the midplane. The second term is a local analog of a logarith- 
mic spiral potential considered in Roberts (1969) and Shu et 
al. (1973). To parametrize the spiral arm strength, we employ 
the dimensionless parameter 



m 

sine 



1$ 



spl 



rI^I 



(7) 



corresponding to the ratio of the maximum force due to the 
spiral potential to the the mean radial gravitational force 
(Roberts 1969). 

The interstellar gas is subject to the net cooling pC = 
n^A[T]-nr, where n = p/{pmn) is the gas number density 
with /i = 1 .27 being the mean molecular weight per particle. 
For the cooling rate of the diffuse ISM, we take the fitting 
formula suggested by Koyama & Inutsuka (2002), 



A(T) = 2x 10"'^exp 



-1.184 X 10= 

r+1000 



-92 



-1-2.8 X 10"^^\/fexp (^^J erg cm^ s"S (8) 

with temperature T in degrees Kelvin. 

For the gas heating function, we consider two different 
cases: (1) a constant heating rate F = Fq = 2.0 x 10"^^ ergs"^ 
and (2) a density-dependent heating rate 



F = Foexp(n/no)^ 



(9) 



with no = 20 cm"^. The first, uniform heating rate corre- 
sponds to the mean heating rate due to the photoelectric effect 
on small grains and PAHs, photodissociations of hydrogen 
molecules, and ionizations by cosmic rays and X-rays (e.g., 
Koyama & Inutsuka 2002). The second, density-modified 
heating rate that increases stiffly with n is to treat the ef- 
fect of star formation feedback in a very simple way. With- 
out such a feedback, high-density clouds produced inside spi- 
ral arms in our self-gravitating models would undergo catas- 
trophic collapse, preventing further computation. In the real 
ISM, gravitational collapse leads to new-born stars which will 
in turn disperse the parent clouds by injecting radiative and 
mechanical energies. Investigating the details of star forma- 
tion and consequent feedback processes is a very active cur- 
rent research area. Most previous work has adopted simplified 
feedback prescriptions that depend on gas-consumption rate, 
star-formation efficiency, type of energy injection, etc., with 
considerable uncertainties in the parameters (e.g., Springel et 
al. 2005; Joung & Mac Low 2006; Shetty & Ostriker 2008; 
Koyama & Ostriker 2009a). More realistic feedback prescrip- 
tions will be considered in a subsequent paper. 

Figure 1 plots the equilibrium cooling curves in the den- 
sity vs. pressure plane. The solid line corresponds to the uni- 
form heating rate, while the dashed curve is for the modi- 
fied heating rate. The dotted lines indicate isotherms. The 
modified heating rate changes the equiUbrium curve dramat- 
ically only for high-density gas, while making a negligible 
difference for low-density material. The equilibrium pres- 
sure has a local maximum Anax/^B = 5.0 x 10^ Kcm"^ at 
icrit.i = 1-0 cm"^ for both heating rates, while attaining a local 
minimum at ncrit.2 = 8.7 cm"^ withP^in/^B = 1-6 x 10^ K cm~^ 
for the constant heating rate, and at ncrit,2 = 6.9 cm"^ with 
Pmm /kn = 1.7 X 10^ Kcm~^ for the modified heating rate. Un- 
der the constant heating rate, the gas temperature along the 
equilibrium curve is a monotonically decreasing function of 
density, although it is insensitive to n at the low-density end 
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log(n) [cm '] 

Fig. 1 . — Thermal equilibrium curves in the density-pressure plane. Solid 
and dashed lines correspond to the uniform heating rate and the density- 
modified heating rate, respectively, while dotted lines indicate isotherms with 
T in Kelvin. 

with n < Merit. 1 • This is not the case under the density-modified 
heating rate where the equilibrium temperature increases with 
density when n > ncrit.3 = (1 /3)'''"'no = 13.9 cm"-' in order to 
model feedback. We thus classify the gas into two compo- 
nents based on its density rather than temperature: rarefied 
component if n < ncrit.i and dense component if n > ncrit,i- 
Note that thermally-unstable gas with «crit,i < « < «crit,2 be- 
longs to the dense phase according to our classification. 

As equation (4) indicates, we explicitly ignore the ef- 
fect of thermal conduction in the present work. Paper I 
found that large translational motions in a finite difference 
scheme give rise to numerical diffusion that tends to sup- 
press the growth of TI, similarly to thermal conduction. The 
amount of numerical conductivity in our models is typi- 
cally /C„ = 10^ ergs"^ cm"' K"' for the background velocity 
vox = 13 kms"\ grid spacing Ax = 2.5 pc, and the perturba- 
tion wavelength A = 20 pc. Inclusion of physical conductivity 
larger than /C„ would resolve the wavelengths of the most un- 
stable Tl. But, this would in turn reduce the timestep greatly, 
making computation essentially unpractical.' We note that by 
neglecting the thermal conduction term in the energy equa- 
tion, some of our results may depend on the numerical reso- 
lution, although the mass fractions appear to be insensitive to 
the resolution (Paper I). 

2.2. Model Parameters & Numerical Methods 

We consider a simulation box located near the Solar neigh- 
borhood at galactocentric radius of /?o = 8 kpc. We adopt 
the galactic rotational velocity of Ro^q = 208 kms"' with a 
flat rotation curve (qo = 1). The corresponding angular ve- 
locity is r^o = 26kms"' kpc"', and orbital period is forb = 
27r/f2o = 2.4 X 10** yr, which we use as the time unit in our 
presentation. For spiral arm parameters, we take pattern speed 
f2p = 0.5^0 > pitch angle sin? = 0.1, and azimuthal wavenum- 
ber m = 2. The corresponding arm-to-arm separation is = 
27r/?osin;/m = 2.5 kpc, which is set equal to the size of the 

' We have run some simulations by including density-dependent thermal 
conductivity IC„ = 10** ergs"' cm"' K"'(l 4-0.05 cm^^/n)"' (Koyama & Os- 
triker 2009a), and confirmed that this level of thermal conduction does not 
make a significant difference in the results. 



simulation box along the x-direction. We fix the spiral arm 
strength to Fq = 5%. 

Our initial gaseous disks, in the absence of the spiral- 
arm perturbations, are taken to be isothermal and in verti- 
cal hydrostatic equilibrium under the Unear stellar gravity 

= —A'KGpt.z (see eq. [6]). The corresponding density dis- 
tribution is a Gaussian profile 

p{z) = poexp (y^j^ , (10) 

with a scale height 

' (47rGp*)i/2 

/ \-l/2 

where cr = 7 km s~' is the isothermal sound speed of the initial 
disks and = 0.056 M0 pc"-' is the stellar density near the 
solar neighborhood (Holmberg & Flynn 2000). We take = 
1 .5hg = 960 pc as the vertical size of the simulation domain 
(i.e., \z\<Lj2). 

Tables 1 and 2 summarize the model parameters and some 
simulation outcomes for models with and without spiral po- 
tential perturbations, respectively. Column (1) labels each 
run. The prefixes NU and NM stand for non-self-gravitating 
models ("N") with the uniform heating rate ("U") and the 
modified heating rate ("M"), respectively, while the prefix 
SM indicates self-gravitating models ("S") with the modified 
heating rate ("M"). As will be discussed below, column (2) 
gives the initial gas surface density Eq. Columns (3) and (4) 
give the mass fractions, = (/ padxdzj J pdxdz) (with a=D 
or R), of dense and rarefied components, respectively. Here, 
the angle brackets ( ) denote a time average over t /forb = 5-8 
for non-self-gravitating models and over f /forb = 8-11 for 
self-gravitating models. Columns (5) and (6) give the scale 
heights, //a= (/ paZ^dxdz/ J Padxi/z)'/^, of the dense and rar- 
efied components, respectively. Column (7) gives the average 

1 /2 

scale height of the whole gas //ave = {foH^ + /rHr) 

We integrate the time-dependent partial differential equa- 
tions (2)-(5) using a modified version of the Athena code 
(Gardiner & Stone 2005, 2008; Stone et al. 2008; Stone & 
Gardiner 2009). Athena employs a single-step, directionally 
unsplit Godunov scheme for magnetohydrodynamics in mul- 
tispatial dimensions. Among the various schemes contained 
in Athena, we take a piecewise linear method for spatial re- 
construction, HLLE Riemann solver to compute the fluxes 
(Harten et al. 1983; Einfeldt et al. 1991), and van Leer al- 
gorithm for directionally unsplit integration (Stone & Gar- 
diner 2009). Since our simulations involve strong shocks for 
the dense medium (with typical Mach numbers ~ 7 - 10), we 
adopt the first order flux correction when the net mass flux out 
of a cell exceeds the initial mass of the cell in order to avoid 
an occurrence of negative density (see, e.g., Lemaster & Stone 
2009). Our models employ a 1024 x 5 12 zones over the simu- 
lation box, corresponding to the resolution of Ax = 2.4 pc and 
Az= 1.9 pc. 

We adopt the shearing-periodic boundary condition at the 
.v'-boundaries (Hawley et al. 1995). In the j-direclion, we use 
the outflow condition for the velocity and the vacuum condi- 
tion for the gravitational potential (e.g., Koyama & Ostriker 
2009a). For the density and pressure at the z-boundaries, 
we linearly extrapolate the logarithmic density, while keep- 
ing temperature fixed, whenever dp/dz < 0. This produces a 
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TABLE 1 

Models Without Spiral Arms {F = 0%) 



Model" 


So (Mo pc-2) 


fo (%) 


fR (%) 


Hd (pc) 


H/f (pc) 


Have (pc) 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


NIT S02 


2 





100 





126 


126 


NU.S05 


5 





100 





119 


119 


NU.SIO 


10 


71 


29 


2 


125 


67 


NU.S20 


20 


86 


14 


4 


127 


48 


NM.S02 


2 





100 





126 


126 


NM.S05 


5 





100 





119 


119 


NM.SIO 


10 


69 


31 


4 


125 


70 


NM.S20 


20 


85 


15 


9 


130 


50 


SM.S02 


2 





100 





121 


121 


SM.S05 


5 





100 





107 


107 


SM.SIO 


10 


82 


18 


4 


100 


43 


SM.S20 


20 


94 


6 


7 


84 


21 



"The prefixes NU refers to the non-self-gravitating models with the uniform heat- 
ing rate, NM for the non-self-gravitating models with the modified heating rate, 
and SM for the self-gravitating models with the modified heating rate. 



TABLE 2 

Models With Spiral Arms (F = 5%) 



Model 


So (M0 pc-2) 


fo (%) 


fR (%) 


Hd (pc) 


Hr (pc) 


H-ive (pc) 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


NU.S02 


2 


11 


89 


25 


129 


122 


NU.S05 


5 


62 


38 


10 


130 


81 


NU.SIO 


10 


81 


19 


7 


130 


57 


NM.S02 


2 


12 


88 


28 


130 


122 


NM.S05 


5 


60 


40 


20 


132 


84 


NM.SIO 


10 


81 


19 


25 


135 


64 


SM.S02 


2 


26 


74 


14 


124 


107 


SM.S05 


5 


91 


9 


21 


121 


42 


SM.S10 


10 


95 


5 


43 


123 


51 



Note. — Model name prefixes are as in Table 1. 



balance between the vertical pressure gradient and the gravita- 
tional source term at the boundaries, similarly to the "conduct- 
ing" boundary in Parrish & Stone (2005). When dp/dz > 0, 
on the other hand, we switch to the continuous boundary con- 
dition for the density and pressure to reduce artificial mass in- 
flow due to the extrapolation. Under our boundary conditions, 
the gas can freely escape from the vertical boundary; we have 
checked that the total mass is nonetheless conserved within 
2% for all models. Because of the very short cooling time, en- 
ergy updates from the net cooling functions are made implic- 
itly based on Newton-Raphson iterations (Piontek & Ostriker 
2004). To solve for the gravitational potential in our simu- 
lation domain, we adopt a method introduced by Koyama & 
Ostriker (2009a) which, by using the fast Fourier transform 
technique, is much more efficient than a hybrid method in- 
volving Green's functions (e.g., Kim et al. 2002) 

3. VERTICAL EQUILIBRIA WITHOUT SPIRAL ARMS 

While our main objective is to investigate the overall dy- 
namics and structure of spiral shocks in verticaUy stratified 
disks under the influence of TI, in this section we focus on the 
quasi-static vertical equilibria with heating and cooling in the 
absence of the spiral arm potential (i.e., F = 0). This allows us 
to study the effect of TI on vertical disk structure. We run one- 
dimensional simulations with physical quantities varying only 
with z. We consider an initially-isothermal disk with Eq = 2, 
5, 10, or 20 M© pc~^, and evolve it subject to TI. 

For disks with large surface density (models with Sq > 
10 Mq pc~^), TI grows rapidly (<C ?orb). transforming the ini- 
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Fig. 2. — Distributions of number density along the vertical direction for 
one-dimensional non-self-gravitating models NU.SIO (solid), NM.SIO (dot- 
ted), and self-gravitating model SM.SIO (dashed) without spiral-arm poten- 
tial perturbations. The midplane densities are n(()) = 32, 15, and 20 cm"' 
for models NU.SIO, NM.SIO, and SM.SIO, respectively. The interface be- 
tween the midplane dense layer and the surrounding rarefied medium oc- 
curs at almost the same density ntrans ~ 0.25-0.35 cm"', corresponding to 
Piims/kB ~ 2000 - 2200 cm-3 K. 

tially constant-temperature gas into thermally bistable phases. 
The cold, dense gas faUs toward the midplane to form a dense 
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Fig. 3. — (a) Mass fraction and (h) scale height of the rarefied compo- 
nent as functions of total surface density So from one-dimensional simula- 
tions without spiral arms. Symbols represent the numerical results for non- 
self-gravitating models NU (asterisks), NM (diamonds), and self-gravitating 
models SM (squares). Solid and dashed curves are the theoretical estimates 
for two-phase equilibrium with and without self-gravity, respectively, for 
which Pi;(0)/k = Ptrans/^B = 2100 K cm"' and cr=1 kms"' are adopted. Ver- 
tical dot-dashed lines mark Sn,in = 2.1 pc"^ and Smax = 6.7 pc"^; 
for Smin < E < Emax, both single-phase and two-phase equilibria are possi- 
ble. 

slab, while the warm, rarefied gas rises up buoyantly. The 
infall is supersonic relative to the dense medium. At early 
time, the dense slab surrounded by the upper rarefied gas 
undergoes vertical expansions and contractions a few times. 
As the kinetic energy dissipates through shocks at the inter- 
faces, the whole configuration evolves toward vertical hy- 
drostatic equilibrium typically within ~ 0.6/orb- Figure 2 
shows density profiles for S 10 models with Eq = 10 pc~^. 
Solid and dotted lines are for non-self-gravitating NU and 
NM models, respectively, while the dashed line is for the 
self-gravitating SM models. The difference between mod- 
els NU.SIO and NM.SIO is not significant since the maxi- 
mum midplane density is not much larger than no = 20 cm"-', 
below which the heating rate is almost density-independent. 
For model SM.SIO, self-gravity compresses the midplane slab 
further at the expense of the rarefied medium at |z| > Hb. 
Nevertheless, the phase transition between dense and rarefied 
components turns out to occur at almost the same density 
itrans ~ 0.25-0.35 cm"^, corresponding to the transition pres- 
sure Ptrans/^B ^ 2000-2200 cm^^ K for all models that are 
unstable to TI. Note that /Vans/^B is above the minimum pres- 
sure for a cold medium with our adopted cooling and heating 
functions, fmi„//tB = 1 600 - 1 700 cm'^ K. 

Once vertical hydrostatic equilibrium is attained, we mea- 
sure the mass fractions fr, and /«, and the scale heights Ho 
and Hr of the dense and rarefied components, respectively; 
these values are listed in Table 1. Figure 3 plots /« and Hr 
as functions of the initial surface density Sq. The results of 
NU and NM models are denoted by asterisks and diamonds, 
respectively, while open squares are for SM models. Mod- 
els with low surface density (So = 2 M© pc"^) do not expe- 
rience TI and thus estabUsh a single-phase equilibrium con- 
sisting only of the rarefied medium. Since the warm gas is 
nearly isothermal at cr « 7kms"^ and self-gravity is weak 
in these models, the equiUbrium density profiles in low S 



models are approximately given by equation (10), with sur- 
face density S = pohgV^ir = p^crI yjTGpl. Since the mid- 
plane pressure /r(0) = cjpo of the rarefied gas cannot ex- 
ceed Pjnax along the thermal equilibrium curve shown in Fig- 
ure 1, the surface density for a single-phase equilibrium with 
only a rarefied component is constrained to be less than 
Smax = Pmax/(2G/9*c|)'''^ = 6.7 M© pc~^. Similarly, the con- 
dition Pr(0) = Pmin yields Smin = 2.1 Mq pc"^ as the minimum 
surface density for a two-phase equilibrium in which dense 
and rarefied components coexist. The two vertical dot-dashed 
lines in Figure 3 mark Smin and Smax- Below Smin. only a rar- 
efied phase is possible, whereas above Smax, both dense and 
rarefied phases must be present. 

For Smin < So < Smax, both siuglc (rarefied) phase and 
two-phase equilibria can be realized. Which type of equi- 
librium emerges depends of course on the initial disk con- 
ditions. In the case of our models with Sq = 5 M© pc~^, the 
initial midplane density and pressure are «(0) = 0.5 cm"^ and 
f (0)/^B = 3770 K cm"-', smaller than than ncrft.i and Pmax/^B- 
Since cooling and heating occur almost isobarically, even the 
densest gas in these models is unable to overcome /\nax to turn 
into the dense component, for this case. 

Figure 3 also shows that for the models that reach a two- 
phase equilibrium, self-gravity reduces the rarefied-gas frac- 
tion in mass as well as its scale height compared to those in 
non-self-gravitating counterparts. Self-gravity also makes the 
density profile of the rarefied component deviate significantly 
from a Gaussian profile. A thin midplane dense slab, contain- 
ing the majority of the gas mass, exerts a uniform gravity on 
the rarefied gas lying above it, providing an additional con- 
fining force. In the Appendix, we describe a simple way to 
estimate //? and Hr as functions of the total gas density, as- 
suming that the rarefied component can be characterized by 
a fixed sound speed cr and that its self-gravity is negligible. 
The resulting theoretical predictions, with cr = 7kms , for 
self-gravitating and non-self-gravitating cases are plotted in 
Figure 3 as dashed and soUd curves, respectively. These are 
overall in good agreement with the numerical results. Small 
discrepancies between the theoretical and numerical values 
of Hr for disks with So = 5 Mq pc"^ arise from the fact that 
the rarefied gas in these models has larger midplane pressure 
than in any other models.^ In view of the thermal equilibrium 
curve shown in Figure 1, this implies that the rarefied medium 
in SOS models is coldest, corresponding to cr ~ 6.3kms"\ 
making the scale height smaller than the theoretical estimate 
based on cs = 7 kms~^ 

4. NON-SELF-GRAVITATING MODELS 

Now we turn to our main theme, nonlinear gas flows with TI 
across spiral arms in a stratified disk. In this section, we study 
overall evolution, structure, and statistical properties such as 
density and random velocity distributions of spiral shocks for 
non-self-gravitating models. Effects of self-gravity will be 
discussed in the next section. 

4.1. Overall Evolution 

We begin by describing evolution of our fiducial models 
NU.SIO and NM.SIO with Sq = 10 Mq pc'^ that employ the 
uniform and density-modified heating rates, respectively. We 

^ For example, the midplane pressure of the rarefied component is 
PR(0)/kB = 3500 and 2500 cm"' K for models NM.S05 and NM.SIO, respec- 
tively. 
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Fig. 4. — Density snapshots for model NU.SIO at ?/foib = 1-50, 2.24, and 2.40. A few instantaneous streamlines are drawn as solid lines in (a). The shock front 
alternates between convex (h) and concave (c) shapes, seen from the upstream direction, due to quasi-peiiodic flapping. Three dense condensations located near 
x/Lx = —0.3, 0, and 0.14 in (b) have moved to x/L, = -0.12, 0.14, and 0.35 in (c), respectively. Colorbar labels log(n/l cm"'). 



slowly turn on the spiral potential amplitude such that it at- 
tains full strength = 5% at f /foib =1-5. Snapshots of volume 
density in logarithmic color scale at early epochs t /foit = 1 -50, 
2.24, and 2.40 are shown in Figures 4 and 5, respectively. 
Figure 6 plots the gas distribution in the n-P/k^ plane for both 
models at f/foib = 2.40. Initially, the disk is in hydrostatic 
equilibrium with a constant sound speed of cr = 7kms"'. 
Since the initial disk is out of thermal equilibrium, it quickly 
evolves into a two-phase equilibrium configuration, as ex- 
plained in §3. As F increases, both the dense gas near the 
midplane and the rarefied gas at high |z| respond to the grow- 
ing spiral potential and form a shock front near the potential 
minimum. 

Since the gas flows at this time are fairly horizontal with- 
out much vertical mixing, as evidenced by the instantaneous 
streamlines shown in Figures Aa and 5a, the shock profile at 
each height is very similar to those of the one-dimensional 
cases studied in Paper II. The shock strength and gas phase in 
the postshock regions depend on the mean density and tem- 
perature at that height. Near the midplane at |z| < Ho (= 7 



and 25 pc for models NU.SIO and NM.SIO, respectively), 
the dense slab is so cold that the shock is very strong with 
a typical Mach number 10, resulting in a far denser 

postshock region. In the high-|z| regions (|z| > 130 pc), on 
the other hand, the gas is warm and has a low mean den- 
sity (< 0.1 cm"^) enough to remain warm even after the 
shock compression. It is the mid-altitude rarefied medium (at 
Hd < \z\ < 130 pc) that is able to achieve a postshock pres- 
sure larger than P^ax and thus undergoes a phase transition 
to the dense component after experiencing isobaric cooling 
(Mufson 1974; Inoue & Inutsuka 2008; Paper II). Since the 
shock is stronger at lower |z| in a stratified disk and since a 
stronger shock tends to move downstream (e.g., Kim & Os- 
triker 2002), the shock front when it first develops is naturally 
curved in the XZ plane. Figures 4a and 5fl show that the shock 
front is concave when seen from the upstream direction, with 
mean slopes of \dxsj,{z)/dz\ ~ 0.83 at |z| < 130 pc and 0.13 at 
|z| > 130 pc, where Xfip(z) is the shock location at z- 

The dense gas produced at the shock at moderate z begins 
falling toward the midplane under the influence of the external 
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Fig. 5. — Density snapshots for model NM.SIO at f/forb = 1-50, 2.2, and 2.40. A few instantaneous streamlines are drawn as solid lines in (a). Compared 
to model NU.SIO, the density-modified heating rate thickens the dense midplane layer and prohibits the formation of dense condensations. Colorbar labels 
log(n/l cm"''). 



gravity as it moves downstream. The reduction of the velocity 
in the direction normal to the concave shock front also helps 
the downward motion of the gas. On the other hand, the dense 
gas near the midplane has a large pressure and thus slightly 
expands vertically after the shock. The vertical expansion is 
more extreme in NM models than in NU models. The falling 
gas collides with the expanding gas, reducing the rising mo- 
tion of the latter. The streamlines shown in Figures 4a and 5fl 
illustrate these motions at early time. 

The rarefied gas which crosses the shock at high |z| also 
falls toward the midplane as it follows galaxy rotation. This 
builds up thermal pressure at low so the flow rebounds to 
high-altitude regions. Since the period of vertical oscillation, 
^ (Gp*)"'/^, is in general not commensurate with the inter- 
val between arm crossings, the streamlines of the rarefied gas 
are not closed. This causes the shock front to sway back and 
forth around its mean position in the direction perpendicular 
to the arm (e.g., Kim & Ostriker 2006; Paper I). During the 
course of the flapping motions, the shock front has a convex 
shape (seen from upstream) when the postshock regions are 



maximally compressed (Fig. 4b), while it becomes concave 
when the gas in the postshock regions is in full vertical ex- 
pansion (Fig. 4c). These flapping motions of the shock front, 
alternating between convex and concave shapes, occur quasi- 
periodically with a period of ^ 0.5forb and have an amplitude 
of Ax/Hr - 1 at |z| = Hr (= 130 pc in model NU.SIO). The 
shock flapping motions are able to tap some of the kinetic en- 
ergy in galaxy rotation to supply random kinetic energy for 
the gas. We will quantify the amplitudes of random gas mo- 
tions driven by flapping in §4.3. 

One of the special features of galactic spiral shocks is that 
gas experiences acceleration after the maximum shock com- 
pression, forming a postshock expansion zone (e.g., Balbus 
1988; Kim & Ostriker 2006; Papers I & II). Any parcel of gas 
becomes gradually less dense as it moves downstream in the 
expansion zone. In model NU.SIO, the shock compression 
and subsequent cooling is so strong that the shocked dense 
gas in the midplane can reach « > 10^ cm"^ (see also model 
SC20 in Paper II). With such a large postshock density, this 
gas can still remain dense, with n > ncrit,2, even after emerg- 
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Fig. 6. — Distribution of gas in the density-pressure plane for models (a) NU.SIO and (b) NM.SIO at f/forb = 2.4, witii grayscale indicating the mass fraction 
in logarithmic scale. The thermal equilibrium curves are the same as in Figure 1. 

ing from the expansion zone located atx/Lx~0-0.3. This 
Tl-stable dense gas travels almost ballistically in the interarm 
region, reenters the arm, and combines with other dense gas to 
produce a few condensations. Figures 4b at time f/forb = 2.24 
shows three large condensations located at x/L^ ^ -0.3, 0, 
and 0.14, which are stretched horizontally due to the expand- 
ing background velocity. The condensations move nearly 
horizontally to appear at x/L^ = -0.12, 0.14, and 0.35 when 
t/t^^ = 2.40 (Fig. 4c). In model NU.SIO, the dense con- 
densations stay in the arm (-0.05 < x/L^ < 0.05) for about 
^ 0.15forb, in the expansion zone (0.05 < x/L^ < 0.25) for 
about ~ 0.30forb. and in the interarm region for the remain- 
der (^ 0.55?orb) of the cycle. The width and residence time of 
each zone are insensitive to the model parameters. 

With the density-modified heating rate, on the other hand, 
the postshock dense gas in model NM.SIO has a moderate 
density (~ 30 - 40 cm"^), so that the postshock expansion is 
able to take it to the thermally unstable regime (ncrit,i <n< 
«crit,2)- Subsequently, the expanding dense gas suffers from TI 
and separates back into dense and rarefied gas in the interarm 
region. The large thermal pressure also prevents the formation 
of dense condensations in this model. 

Figure 7 plots the temporal evolution of the mass fractions 
of dense phase {solid lines) and rarefied phase {dashed lines), 
respectively, for models NU.SIO {thick lines) and NM.SIO 
{thin lines). At early time, fo increases rapidly as the gas 
cools and collapses toward the midplane to form a dense 
layer that bounces appreciably at ///orb ~ 0.1. The mass frac- 
tions flatten at ///orb ~ 0.6 when vertical hydrostatic equilib- 
rium is established, well before the effect of the spiral po- 
tential becomes substantial. As the spiral potential attains 
its full strength at ///orb = 1-5, fo increases slightly due to 
the phase transition of the rarefied to dense phases occurring 
at the shock fronts. Although the flows are fully nonlinear 
with strong unsteady motions and phase transitions, there is 
no noticeable secular variation in the mass fractions, which 
remains at fo ~ 0.8 after ///orb = 5; the associated tempo- 
ral fluctuation amplitudes are about 6-9% relative to the 
mean values. We thus conclude that in a statistical sense, 
the spiral shocks in our models have reached a quasi- steady 
state at ///orb > 5. Compared with models without spiral 



Fig. 7. — Mass fractions of the dense (fo) and rarefied (/r) components 
as functions of time for models NU.SIO {thick) and NM.SIO {thin). Initially, 
fo increases rapidly as the gas cools due to Tl and flattens to fo ~ 0.7 at 
t/toih ~ 0.6 when hydrostatic equilibrium is attained before the effect of the 
spiral potential becomes significant. The presence of the spiral arm at full 
strength increases this to a saturate value of fo ~ 0.8 at f/f„,|, ^ 5. The mass 
fractions of model NU (unmodified heating) and NM (modified heating) are 
quite similar. 

arms discussed in §3, the shock compression and associ- 
ated phase transitions decrease the rarefied gas fraction by 
46% for SIO models. In fact, all of the non-self-gravitating 
models with spiral arms have comparable total surface den- 
sity of rarefied gas, ~ 1.9 M© pc"^, lower than the value 
- PtTms/{CR\/2Gp») « 2.8 M© pc~^ that would be pre- 
dicted using a uniform surface density. Note that both NU 
and NM models have almost the same dense and rarefied mass 
fractions since the modified heating rate does not affect the 
rarefied medium much. 

The evolution of S02 and SOS models is quaUtatively simi- 
lar to that of SIO models in that phase transitions occur at the 
shock and in the postshock expansion zone, although the for- 
mer with low postshock density do not produce much dense 
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Fig. 8. — Density distributions averaged over ;/<orb = 5-8 of our non-self-gravitating models. S05 and SIO models contain midplane dense gas in both arm 
and interarm regions, while S()2 models have the dense phase only in the arm regions. Compared to NU models, the arm regions in NM models are broader and 
thicker. Colorbar labels log(n/l cm"'). 



gas even under the uniform heating rate. When the spiral po- 
tential is absent, the equilibrium disks of these models con- 
sist entirely of the rarefied gas with the midplane pressure 
P(0) /ks - 3500-4000 cm'^ K for 805 models and P(0) /k^ - 
1500- 2000 cm"^^ K for S02 models. But, the shock com- 
pression increases the pressure by about a factor of 3, cor- 
responding to a typical Mach number w vo.v/cr ^ 2 for the 
rarefied medium,"* making the midplane postshock pressure 
larger than Pmax- As a result, the dense medium in S05 mod- 
els comprises about 60% of the total mass, undergoing Tl in 
the postshock expansion zone. In S02 models, the postshock 
pressure barely exceeds Pmax, so that the shocked dense gas, 
comprising about 10% of the total, easily disperses to return 
to the rarefied gas in the interarm region. Flapping motions 
of spiral shocks are correspondingly weaker in these models, 
with Ax/Hr ^ 0.7 and 0.3 for S05 and S02 models, respec- 
tively. 

4.2. Time-averaged Shock Structure 

To visualize spiral shock structure in each model, we con- 
struct a time-averaged distribution of number density (n). 
Here, the angle brackets ( ) denote a time average over t / to^b = 
5-8. Figure 8 displays (n) for all the non-self-gravitating 

' For adiabatic shocks, the pressure jump condition is P2/P1 = 1 + A4"(l — 
l/s), where the subscripts 1 and 2 denote preshock and postshock values, 
respectively, andi= [(■y+l)M^]/[2'y+(j-l)M^] is the density shock jump 
factor 



models in logarithmic color scale. It is apparent that S05 and 
SIO models possess a thin dense layer everywhere near the 
midplane, while the dense gas is found only inside the arm 
regions in S02 models. The shock compression factors in the 
time-averaged density profiles are ^7-10, which is larger 
than the adiabatic shock jump due to enhanced radiative cool- 
ing in the shocked gas (cf., Mufson 1974; Inoue & Inutsuka 
2008, Paper II). The shock transition layer in S05 and SIO 
models is relatively broad because of rather strong flapping 
motions of the shocks, while S02 models exhibit relatively 
sharp discontinuities. Compared to NU models, arms in NM 
models have larger pressure and are more expanded vertically, 
similar to "hydraulic jumps" that occur when the equation of 
state is stiffer than isothermal (e.g., Martos & Cox 1998). Ta- 
ble 2 lists the time-averaged values of the mass fractions and 
scale heights of dense and rarefied components, as well as the 
overall average scale height. 

Figure 9 plots the mass-weighted probability distribution 
functions (PDFs), averaged over f/forb = 5-8, of gas den- 
sity and temperature for models NU.SIO and NM.SIO. The 
PDFs are in general bimodal, as is expected for a bistable 
cooling function. For model NU.SIO, the dense and rar- 
efied peaks are centered at {n,T) ^ (200 cm""* , 30 K) and ^ 
(0.2 cm"-', 7100 K), respectively, mostly distributed near the 
thermal equilibrium curves. Only a small fraction of the gas 
is in the thermally-unstable regime. The dense portion of the 
PDF in model NU.SIO is well fitted by a lognormal distribu- 
tion {thin solid line) with a standard deviation of A(lnn) = 1 .2, 
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Fig. 9. — Mass-weighted density (top) and temperature (bottom) prob- 
ability density functions, averaged over t/t^ib = 5-8, in models NU.SIO 
(solid) and NM.SIO (dotted). While the broad dense peak centered at 
(n.T) ^ (200 cm"' , 30 K) in model NU.SIO is compressed and shifted 
to (20cm-3.180K) in model NM.SIO. the rarefied peak at (n.T) ~ 
(0.2 cm"', 7100 K) is unchanged. The thin line in the top panel is a log- 
normal fit to the dense peak in model NU.SIO, with a standard deviation of 
A(lnn) = 1.2. 

which is one of the characteristics of near-isothermal turbu- 
lence (cf., Wada & Norman 2007; Wada 2008). With en- 
hanced heating, on the other hand, model NM.SIO shows a 
sharp density cutoff in the density PDF at n ^ 50 cm""* and has 
a dense peak shifted to n ~ 20 cm"^. Because of the stiff equa- 
tion of state, the dense gas in model NM.SIO is not as cold as 
in model NM.SIO. This not only thickens the midplane dense 
slab, but also sets an upper limit on the gas density, which in 
turn prevents the formation of dense condensations. In model 
NM.SIO, all the postshock dense gas becomes thermally un- 
stable in the expansion zone and separates into dense and rar- 
efied gas in the interarm region. 

Paper II showed that for one-dimensional models, the den- 
sity profiles of multiphase spiral shocks are more symmetric 
and have a wider arms than isothermal counterparts. This is 
because the strength of spiral shocks in the multiphase models 
fluctuates depending on whether the incoming gas is warm or 
cold, resulting in slight oscillations of the shock fronts in the 
direction perpendicular to the arm. In addition, spiral shocks 
in the XZ plane undergo flapping motions, which can further 
widen the arms. To see this, we plot in Figure 10 the time- 
averaged surface density profiles = J_^{p)dz after tak- 
ing a boxcar average with window of 8 pc. The solid and 
dotted lines are for NU and NM models, respectively. Shown 
also as dashed lines are the density profiles n(x)/no from one- 
dimensional simulations (i.e. without vertical stratification) 
under the uniform heating rate; the initial number density no 
of the one-dimensional counterpart was chosen equal to the 
density-weighted mean density Have = 2]o/(27r'''^/imH//ave), 
with the average disk thickness i/ave listed in column (7) of 
Table 2. For SOS and SIO models for which the shock flap- 
ping motions are appreciable, the arms are considerably wider 
and less centrally peaked than in the one-dimensional mod- 
els. Due to the flapping motions, dense condensates formed in 
tiie NU.SIO model oscillate slightly in the jc-direction when 
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Fig. 10. — Time-averaged profiles of surface density from our two- 
dimensional simulations under the uniform heating rate (solid) and the 
density-modified heating rate (dotted) for models with (a) Eq = 10 pc"^, 
(b) 5 Mq pc"^, and (c) 2 Mq pc"^. Dashed lines give the results of one- 
dimensional models (which do not have shock flapping) with the uniform 
heating rate. Stronger shock flapping motions in two-dimensional models 
make the arms broader and less peaked compared with the one-dimensional 
counterparts. 
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Fig. 11. — Veilical scale heights of the time-averaged density distribu- 
tions in (a) NU and (b) NM models as functions of x. The simulation results 
(solid) are overall in good agreement with the theoretical estimates (dotted) 
for effective vertical hydrostatic equilibrium. 

they pass through the shock, resulting in broader arms than 
in the NM.SIO model. For S02 models, the one-dimensional 
shock consists only of the warm rarefied gas, while the two- 
dimensional shocks contain a small amount 10%) of dense 
gas due to an additional compression in the z-direction; the 
difference in profiles is therefore slight. 
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In studies of galactic disk structure, it has been the custom- 
ary to assume effective hydrostatic equilibrium in the vertical 
direction. Using numerical simulations without spiral arms, 
Koyama & Ostriker (2009b) explicitly demonstrated that tur- 
bulent, multiphase disks are in effective hydrostatic equilib- 
rium, provided that the turbulent pressure arising from ran- 
dom gas motions is taken into account. When a spiral poten- 
tial is present, the gas surface density and velocity dispersions 
depend upon the distance x from the minimum of the spiral 
potential. It is interesting to study whether "local" effective 
hydrostatic equilibrium is still established at each x. 

From the time-averaged density distribution, we measure 
the density-weighted vertical scale height H{x), sound speed 
Cs{x) , and vertical velocity dispersion Su^ix) via 



rarefied 



dense 



j {p)dz 



and (5m (x) 



' I{p)dz' 
{vzmdz 



J{p)dz 



(12) 



In the absence of self-gravity, the "estimated" vertical scale 
height is given by 



c1 + 5ui 



(13) 



for effective hydrostatic equilibrium (Koyama & Ostriker 
2009b). Figure 1 1 plots H{x) {solid lines) and i/est(-«) {dotted 
lines) for NU and NM models as functions of x. The mea- 
sured vertical scale height is overall in excellent agreement 
with the estimated value at all horizontal locations. For all 
models, Cj is about 5-7 times larger than 5u^. This implies that 
the disks with spiral arms, in time-averaged sense, are effec- 
tively in vertical hydrostatic equilibrium, with the main sup- 
port provided by thermal pressure (in these models without 
stellar feedback). Unsteady motions associated with shock 
flapping and movements of dense condensations are mostly 
horizontal, affecting the vertical force balance relatively little. 
A small mismatch between H and //est near x/L, = -0. 1 in S02 
models arises from the fact that shocks in these models exhibit 
weak flapping motions and retain sharp discontinuities in the 
time-averaged configurations. In this case, the d{pVxV:^)/dx 
term in the momentum equation has a non-negligible contri- 
bution to the vertical force balance, which was ignored in the 
derivation of equation (13). 

4.3. Velocity Dispersions 

Paper 1 showed that two-dimensional (isothermal) spiral 
shocks exhibit strong flapping motions in the XZ plane and 
are able to generate a sonic level of random gas motions in 
the arm regions. On the other hand. Paper II showed that in 
one-dimensional spiral shocks with TI, random gas motions 
amount to only ^ 1 -2 kms"'. In this subsection, we quan- 
tify the level of random gas motions driven by shock flapping 
motions and TI in our two-dimensional models. 

The velocity field of gas moving across spiral arms is 
a combination of several different components including 
streaming motions, oscillations of the shock fronts them- 
selves, and random motions. Since streaming velocities that 
are ordered and vary with x are much larger than the true ran- 
dom motions of the gas, it is important to subtract the former 
from the total velocity as cleanly as possible. For this purpose, 
we first construct time-averaged templates of the velocity field 
(v,) (with i = X, y, or z) for the dense and rarefied components 
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Fig. 12. — Temporal changes of the density-weighted velocity dispersions 
Svx, Svy, and Sv, of the rarefied (left) and dense (right) components in models 
NU.SIO (solid) and NM.SIO (dashed). The large-ampUtude fluctuations of 
the velocity dispersions, with periods of ~ O.Sforb, are due to incomplete 
subtraction of the arm streaming motions associated with the shock flapping. 

separately. We then calculate the density-weighted velocity 
dispersions using 



Svi{t) = 



J p[vi-{vi)fdxdz 



J pdxdz 



1/2 



(14) 



Figure 12 plots (5v,(f) for the dense and rarefied components in 
models NU.SIO and NM.SIO as solid and dashed curves, re- 
spectively, for a time span of f/^orb = 5-8. Figure 13 plots 
the mean values (Jv,) along with their standard deviations 
A6vi = {{Svj) - (5v,)^)'/^ for all the non-self-gravitating mod- 
els; the values of (Svi) and ASvi are listed in Table 3. 

Figure 12 shows that the density- weigh ted velocity disper- 
sions for the dense component exhibit large-amplitude fluc- 
tuations, with periods roughly of ^ O.Sforb- The standard de- 
viations of the fluctuations are AJv, ~ (0.2-0.5)((5v,) for the 
dense component; deviations are only A^v, ~ (0.1-0.2)((5v,) 
for the rarefied phase. These variations of Sv^ and Svy are 
caused mostly by oscillations of the shock front relative to the 
mean position. With large spatial variations of streaming ve- 
locities across the arm, the small offset of the shock position 
as well as the instantaneous locations of dense condensates 
result in large values of Ai5v,. We thus regard the local min- 
ima of 6vi, approximately equal to cr, = (Svj) - A(5v,, as the 
upper limits to the level of random gas motions. 

Figure 13 shows that for NU models, (Svx) and {5vy) in- 
crease with Eq. This is mainly because the shock compression 
and associated phase transition are stronger with larger Sq, 
leading to stronger flapping motions. Nevertheless, cr^ CTj, ~ 
2-3 kms~^ for both dense and rarefied components, insensi- 
tive to Sq. This indicates that the portion of kinetic energy 
in the shock flapping motions that goes into random gas mo- 
tions is quite limited. The remaining portion is simply associ- 
ated with the horizontal shock oscillations near the midplane. 
Since the shock flapping motions at low |z| are mostly hor- 
izontal, the random vertical motions of the dense gas in NU 
models are forced predominantly by the impact of rarefied gas 
arriving from high altitudes. The fact that the rarefied gas has 
(Tj ~ 1.7 kms~', almost independent of Sqj suggests the flap- 
ping motions drives more-or-less constant vertical motions at 
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TABLE 3 

Induced Random Velocity Dispersions 







Dense Component 




Rarefied Component 


Model 


So 
















(Mq pc-2) 


(kms-') 


(km s-') 


(km s"') 


(kms-') 


(kms-') 


(km s-l) 


NU.S02 


2 


3.04± 1.53 


3.28 ± 1.23 


1.44 ±0.34 


2.62 ±0.58 


2.25 ±0.29 


1.61 ±0.17 


NU.S05 


5 


3.12±0.87 


3.20 ±0.92 


0.61 ±0.14 


3.04 ±0.34 


2.75 ±0.33 


1.64 ±0.14 


NU.SIO 


10 


4.26 ±1.75 


3.60± 1.89 


0.41 ±0.08 


3.47 ±0.81 


2.94 ±0.41 


1.72±0.14 


NM.S02 


2 


3.07 ±1.24 


3.24 ± 1.20 


1.52 ±0.43 


2.79 ±0.40 


2.45 ±0.18 


1.82±0.11 


NM.SOS 


5 


2.99 ±0.66 


3.20 ±0.81 


1.03 ±0.29 


3.17±0.59 


2.69 ±0.46 


1.74±0.11 


NM.SIO 


10 


3.18±0.62 


3.51 ±1.08 


1.11 ±0.40 


3.07 ±0.42 


2.63 ±0.33 


1.69±0.12 


SM.S02 


2 


4.47 ±1.33 


3.97 ±1.49 


0.72 ±0.18 


3.04 ±0.50 


2.49 ±0.19 


1.85±0.14 


SM.S05 


5 


6.36 ±1.79 


6.95 ±2. 11 


0.75 ±0.10 


6.07 ±1.44 


5.50 ±1.45 


2.20 ±0.14 


SM.SIO 


10 


10.52 ±5.55 


8.41 ±3.91 


3.60 ±1.71 


10.60 ±2.94 


7.92 ±1.78 


4.96 ±1.26 



Note. — Model name prefixes are as in Table 1. 
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Fig. 13. — Mean values (symbols) and standard deviations (errorhars) 
of the density-weighted velocity dispersions, averaged over ?/?oib = 5-8, of 
the rarefied (left) and dense (right) components in all non-self-gravitating 
models under the uniform heating rate (top) and the density-modified heating 
rate (bottom). Allowing for the incomplete subtraction of the arm streaming 
motions, the random gas motions are ax ~ o"y ~ 2 — 3 kms"' for both dense 
and rarefied components, and ct^ ~ 1.7 km s"' for the rarefied component, 
insensitive to Sq, in both NU and NM models, while (Tj 
dense component in NU models. 

high |z|; this is because the total mass of rarefied gas is al- 
most the same in all models, equivalent to a surface density 
of 1.9 M0 pc~^. Since the fraction of the rarefied component 
decreases with increasing So (see Table 2), the ratio of ver- 
tical kinetic energy in the rarefied gas to the mass of dense 
gas decreases with Sq. This causes of the dense medium 
to decrease with increasing Eq, roughly as cx Sq"^ 
models. For NM models, the dense gas in the immediate post- 
shock region is overpressured due to the strong heating and 
thus expands vertically, enhancing compared to those in 
NU models. 

5. SELF-GRAVrrATING MODELS 

We now explore the formation of self-gravitating clumps 
and their internal properties. For NU models, we find that 
the inclusion of self-gravity always results in catastrophic col- 
lapses of self-gravitating clouds that form in the postshock re- 
gion, preventing us from continuing simulations further. For 
this reason, we present the results of self-gravitating mod- 



els only with the density-modified heating rate ("SM" mod- 
els). Instead of running self-gravitating models from f = 0, we 
make use of the "saturated-state" data sets of NM models at 
t/torb ^ 4.8 and restart them by slowly turning on the gaseous 
self-gravity over a time interval of 1 .5forb- 

Neglecting the effect of the rarefied medium, the gravita- 
tional susceptibility of a midplane dense layer in NM models 
can be measured by the average Toomre stability parameter 

-I 



Qd 



1 KCd 



0.27 



So 



(15) 



where c/j = 1 kms"' is the mean sound speed of the dense gas. 
With the inclusion of self-gravity, the dense layer in model 
SM.SIO has Qd = 0.34 and thus is quite unstable, initiating 
the collapse of high-density regions. Due to the stiff equation 
of state, however, the collapsing clouds soon reach an equi- 
hbrium state with only a moderate central density 30 cm"^, 
which is only 1 .5 times larger than the average density of the 
dense gas in model NM.SIO. As these clouds follow galaxy 
rotation, they merge together in regions of converging stream- 
ing velocities, eventually resulting in two big condensations. 
Figure 14 shows the density snapshots in logarithmic color 
scale of model SM.SIO at t/toA, = 7.32, 7.48, and 7.69. Two 
clouds are widely separated from each other during traver- 
sal of the interarm region (Fig. 14a). After one cloud en- 
ters the spiral shock, it loses most of its x-momentum, and 
the second cloud can then collide with it when it enters the 
shock at high speed (Fig. \Ab). Since the two clouds are on 
their own epicycUc orbits before the collision, however, they 
retain quite different Vy even after the collision, preventing 
them from merging into a single entity. Due to the Coriolis 
force, the two clouds subsequently have different v^, so that 
the merged entity elongates in the expansion zone (Fig. 14c), 
and separates back into two pieces in the interarm region. 

The dense layer in model SM.S05 has Qd = 0.9 and is thus 
marginally unstable. With the aid of shock compression, the 
dense gas in the postshock region collapses and eventually 
forms two self-gravitating clumps. In this model, the collision 
of these clumps at the shock and subsequent break up in the 
expansion zone is similar to model SM.SIO. With Qd = 12, 
on the other hand, the dense layer in model SM.S02 is gravi- 
tationally stable and does not form dense clumps. Compared 
with model NM.S02, self-gravity in model SM.S02 increases 
the fraction of the dense phase by about a factor of 2, which 
in turn decreases its vertical velocity dispersion by a similar 
factor. 

The presence of self-gravity leads to stronger shock flap- 
ping motions than in the NU models, increasing the veloc- 
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Fig. 14. — Density snapshots of self-gravitating model SM.SIO at ?/?orb = 7-32, 7.48, and 7.69 in logarithmic color scale. Two clouds are separate from each 
other in the interarm regions (a), temporarily merge together at the shock (b), and then break up back into two pieces in the postshock expansion zone (c). The 
boxes surrounding the clouds A and B in (a) are enlarged in Figure 15 to show the interval velocity structure. Colorbar labels log(«/l cm"'). 



ity dispersion. In model SM.S02, the in-plane velocity dis- 
persions of the dense component increases by a factor of 
^ 1.2-1.5 compared to model NU.S02; for low surface den- 
sity, the self-gravity is not sufficiently strong to have a major 
effect. For the SM.S05 and SM.SIO models, however, the 
stronger self-gravity make a larger difference to shock flap- 
ping, in turn driving larger velocity dispersions. Correcting 
for streaming, the velocity dispersions of both dense and rar- 
efied phases in model SM.S05 reach ^ o-y ^ 4-5 kms"', 
about twice larger than those in model NM.S05. In model 
SM.SIO, the dense gas velocity dispersions are ^ cr,. ^ 
4-5 kms"', while the rarefied gas has the velocity dispersions 
up to ^ 7 kms"'. Note that these in-plane velocity dispersions 
in multi-phase, self-gravitating models are similar to those in 
the isothermal self-gravitating models with F = 5% studied in 
Paper I. 

Since the self-gravitating clumps produced in models 
SM.S05 and SM.SIO move almost baflistically, the position 
of the largest surface density at a given time does not always 
correspond to the minimum of the spiral potential. In fact, the 



clumps are near the potential minimum (\x\/Lx <Q.l) only for 
/ forb = 0.35, making the definition of the arm regions rather 
ambiguous. In addition, due to accretion onto the clumps, the 
rarefied medium in these models amounts to less than 10% of 
the total mass, much smaller than the observed mass fraction 
of the warm gas near the solar neighborhood (Heiles 2001; 
Heiles & Troland 2003). For these reasons, these clumps are 
unlike real self-gravitating clouds in spiral galaxies. Never- 
theless, we believe these model clouds may provide clues to 
the internal properties and virial balance of real interstellar 
clouds, in that they represent a limiting situation in which in- 
ternal turbulent feedback from star formation is absent. 

Keeping in mind the caveats mentioned above, we pro- 
ceed as follows to calculate the cloud properties. To define 
the boundary of the dense clumps, we choose the thresh- 
old density nth = 21 cm""*, corresponding to P^ax in the ther- 
mal equilibrium curve. Using a CLUMPFIND algorithm 
(e.g., Williams et al. 1994), we find the interior of each 
cloud with n > nth. We then measure the mean density 
Pel = p, the averaged sound speed Cci = {P/ p)^/^ and the mean 
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TABLE 4 

A\ LRAf;ii PROPIiRI ILS ()!■ CHJMl'S PRODDCLD IN SLLI -GR.W I I ATING MODLLS 



Model tci (Tel Rc\ «d Mc\ a 

(km^-l) (kmi-') (pc) (cm'^) (10*' Mq) 

SM.S05 3.69 ±0.67 0.96 ± 0.24 44.9 ±7.6 27.4 ±0.8 0.36 ±0.15 2.34 ±0.25 

SM.S10 5.06±0.63 1.53±0.61 60.7±6.0 29.8±1.6 0.91 ±0.29 2.31 ±0.32 



one-dimensional internal velocity dispersion CTcI = ^,(*^'? ~ 
Vi^)'/^/3'/^ of each cloud, where the overlines denote the 
mass- weighted average. We then count the total number of 
pixels A'ci on the XZ plane occupied by each cloud, and calcu- 
late the cloud size Ra = (A^ciAxAz/Tr)'/^. Assuming a spher- 
ical shape, we calculate the total mass M^i = 47rpci^ci/3 and 
the virial ratio of each cloud via 

o.Mlfl^l. ,,6, 
GMa 

which is the ratio of the total kinetic energy to the gravita- 
tional potential energy for a uniform spherical cloud (e.g., 
Bertoldi & McKee 1992; McKee & Ostriker 2007); note that 
central concentration would decrease a. Figure 15 gives an 
example of the shapes and internal velocities of two clumps 
shown in Figure 14a. Note that the velocity vectors are dis- 
tributed quite randomly, indicative of (subsonic) turbulent 
motions within the clouds. 

Table 4 summarizes the average properties of the clouds 
that form in models SM.S05 and SM.SIO. The typical size 
and mass of the clouds are found to be Rc\ ^ 45-60 pc and 
Mci ~ (4-9) X 10^ Mq, respectively, with the clouds in model 
SM.S05 somewhat smaller and less massive than in model 
SM.SIO. Overall, a ^ 2 for all the clouds, suggesting that 
they are (marginally) gravitationally bound. The mean sound 
speed inside the clouds is Cd = 3.7-5.1 kms"', ~ 3-4 times 
larger than the one-dimensional internal velocity dispersions 
fTci = 0.9- 1.8 km s~'. This indicates that the major support 
against self-gravity comes predominantly from the thermal 
energy, a consequence of the density-modified heating rate 
we adopt. The relatively low value of the internal turbulent 
velocity dispersion suggests that the interaction of a dense, 
gravitationaUy-bound cloud with its surroundings can drive 
only a modest level of internal turbulence. 

6. SUMMARY & DISCUSSION 



While stellar spiral arms in disk galaxies provide smoothly 
varying low-amplitude gravitational potential perturbations, 
the response of the interstellar gas to them is quite dramatic. 
Spiral shocks compress the ISM and the high post-shock den- 
sities may trigger growth of arm substructures and star for- 
mation. In addition, radiative cooling and heating of the gas 
makes the ISM inherently inhomogeneous, producing two 
phases that differ in density and temperature by about two 
orders of magnitude. Moreover, the vertical stellar gravity 
tends to produce stratification of the cold and warm gas due 
to differential buoyancy; this stratification can be modified by 
vertical turbulence, however. Interactions of these processes 
may significantly affect the gas flows and shock structures, 
compared with results from our previous work (and that of 
other groups), which employed an isothermal approximation 
(Paper I) or neglected the vertical degree of freedom (Pa- 
per 11). In this paper, we have conducted nonlinear hydro- 
dynamic simulations in a two-dimensional slice perpendicu- 
lar to a local segment of a spiral arm that is tightly wound 
(with a pitch angle sin; = 0.1) and rotates rigidly (with a pat- 
tern speed ilp/^o = 0-5). To handle the Coriolis force arising 
from the galaxy rotation self-consistently, we allow for gas 
motions parallel to the arm (i.e. perpendicular to the domain 
of the simulation). We consider two different forms of gas 
heating; the usual constant heating rate (for NU models) and 
the density-modified heating rate (for NM and SM models), 
which mimics the effect of star-formation feedback in a very 
simple way, to limit runaway collapse. We start from initially 
isothermal disks that are in vertical hydrostatic equilibrium 
but out of thermal equilibrium. We slowly turn on the am- 
plitude of the spiral arm potential such that it attains a fuU 
strength at 1.5 orbital times. Magnetic fields are neglected in 
the present work. 

Our main results and their astronomical implications are as 
follows. 

1. Two-phase disk equilibria without spiral arms. — 
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In the absence of spiral-arm potential perturbations (and 
other sources of turbulence), the vertical structure of equilib- 
rium disks depends on the disk surface density So- When 
So > S„>ax = Pmax/(2Gp*4)l/2 ^ 6.7 Mq ^c'^ (for So- 
lar neighborhood conditions), the disk experiences TI and 
evolves toward an equilibrium configuration with a thin slab 
of dense gas (n > 1 cm"^) near the midplane sandwiched 
between layers of rarefied gas (with n < I cm~^). Here, 
Pmax/A^B = 5000 cm~^ K is the maximum pressure allowed for 
the thermally-stable rarefied gas with our adopted heating 
and cooling functions, cr « 7 kms~' is its density-weighted 
sound speed, and = 0.056 M© pc"-' is the stellar density 
near the Solar neighborhood. In our models, the transi- 
tion between the dense and rarefied phases occurs approxi- 
mately at -Ptrans/^B ~ 2100 cm"-' K, insensitive to So, as long 
as the disk is unstable to TI. Without self-gravity, the verti- 
cal distribution of the rarefied gas in equilibrium is well fit- 
ted by a Gaussian profile whose surface density is fixed to 
T,R = Sng = Ptrans/(2Gp*4)i/2 2.8 Mg pc-l When self- 
gravity is included, the gravity from the midplane dense layer 
compresses the overlying rarefied component further, forcing 
overpressured rarefied gas to transform to the dense phase. 
The resulting surface density of the rarefied gas is reduced to 
T,R = /"Sng, with the reduction factor J" defined by equation 
(A5). When S,, < S^n = P„,„/(2Gp,4)'/2 ^ 2.1 Mq pc'^, 
the disk has too low a pressure to produce the dense com- 
ponent; equilibrium consists only of the rarefied gas. When 
Smin < So < Smax. either a single rarefied disk or a two-phase 
disk is possible, depending on the initial conditions. Our S02 
and S05 models with So = 2 and 5 Mq pc"^, respectively, that 
start from a warm isothermal configurations all end up with a 
single-component rarefied disk, in the absence of spiral per- 
turbations. 

2. Shock flapping motions in vertically stratified disks. — 
Spiral-arm potential perturbations lead to spiral shocks in the 
gas, which are vertically curved and non-stationary, show- 
ing strong flapping motions perpendicular to the arms. Sim- 
ilarly to the one-dimensional cases studied in Paper II, the 
shock compression and postshock expansion in two dimen- 
sions allow phase transitions, but only if the gas density at 
the shock and/or the postshock expansion zone reaches the 
thermally-unstable range (1 cm"^ ^ 7-9 cm"^). In model 
NU.SIO with So = 10 M© pc"^ and the uniform heating rate, 
the shocked dense gas has large enough density that the post- 
shock expansion is unable to return it to the thermally un- 
stable regime. As a consequence, there is a large amount of 
interarm dense gas entering the shock in this model, which 
colhdes with other dense gas in the arm, producing dense con- 
densations. In other models with lower surface density or the 
density-modified heating rate, the shocked gas re-expands and 
becomes thermally unstable, returning to either the dense or 
the rarefied phase in the interarm region. 

The shock flapping motions in our models arise due to the 
fact that the arm crossing time of gas is incommensurate with 
the vertical oscillation period, so that steady flows are not 
possible. Seen from the upstream side, the shock is convex 
when the postshock regions are maximally compressed, and 
concave when the postshock vertical expansion is strongest. 
These non-steady motions of shock fronts are conunonly seen 
in numerical simulations with sufficient resolution (e.g., Mar- 
tos & Cox 1998; Gomez & Cox 2002, 2004; Kim & Ostriker 
2006; Paper I; Wada 2008), although simulations with low 
resolution (Tubbs 1980) or particles (e.g., Dobbs et al. 2008) 



do not capture the flapping motions clearly. 

It is interesting to note that radio continuum images of the 
5-kpc arm (or the Scutum arm) toward the galactic longi- 
tude / = 30° - 32° in the Milky Way shows a bow shock fea- 
ture in the warm ionized medium with temperature ~ lO'* K 
(Sofue 1985), similarly to a convex shock front seen in Fig- 
ure 4b. The radio emission from the bow shock is presum- 
ably thermal radiations from ionized gas, with emission mea- 
sure of ~ 7,000 pc cm"^. The curvature of the observed bow 
shock, as measured by the longitudinal offset Al of the shock 
at latitude b relative to the shock front at the midplane, is 
Al/h ^ 0.5 for/? = 0.5°. This value is about a half of the maxi- 
mum curvature of the shock front |xsp(//r)-Xsp(0)|///r ^ 0.85 
in our SIO models, where Xsp(z) denotes the shock position at 
height z. This strongly suggests that the bow shock associ- 
ated with the 5-kpc arm is most likely a cross section of a 
galactic spiral shock that is undergoing flapping motions. Ve- 
locity information is needed to determine whether the 5-kpc 
arm regions are currently being compressed or expanding in 
the course of the flapping motions. 

In this paper, for consistency with our local approximation 
we have considered tightly wound arms with a very small 
pitch angle / ^ 5.8°, and a pattern speed ilp/^o = 0.5. Ob- 
served spiral arms of external galaxies are often more loosely 
wound with / - 10° -30° (e.g., Seigar et al. 2008) and span 
a wide range of galactocentric radii with differing flp/Qo. 
For fixed F, a larger arm pitch angle would imply a larger 
streaming velocity vqx perpendicular to the arm (see eq. 1). 
Spiral shocks would then become correspondingly stronger 
and shifted farther downstream (e.g., Roberts 1969; Shu et 
al. 1973; Kim & Ostriker 2002), exhibiting larger amplitude 
flapping motions (Paper 1). On the other hand, | vox | is increas- 
ingly small as approaches fio. Consequendy, the spiral 
shock as weU as associated flapping motions would become 
weaker as one approaches corotation, where the gas would 
simply concentrate near the potential minimum, without in- 
volving a shock. 

3. Time-averaged shock structure. — Within a few orbital 
times after the development of spiral shocks, gas flows reach a 
quasi-steady state in the sense that the mass fractions of dense 
and rarefied gas do not change appreciably with time. For 
models with So = 10 Mq pc~^, the quasi-steady mass fraction 
of the rarefied gas is ^ 19%, which can be compared to 
/r ^ 30% when the spiral potential is absent. Despite the 
shock flapping motions, most of the gas is found close to 
thermal equilibrium, with a small fraction thermally unsta- 
ble. The density and temperature PDFs are thus bimodal. For 
model NU.SIO, the dense and rarefied peaks are located at 
{nj) - (200cm-3,30K) and - (0.2 cm-^7100 K), respec- 
tively, and the dense part of the density PDF is described 
by a lognormal distribution. The time-averaged structure can 
be well represented by local vertical hydrostatic equilibrium, 
supported mainly by the thermal pressure rather than gas ran- 
dom motions. This indicates that the vertical hydrostatic bal- 
ance is a reasonable approximation even in the presence of 
spiral shocks. The profiles of surface density perpendicular 
to the arm are more-or-less symmetric with a shock compres- 
sion factor of ~ 7-10, and have broad arm regions whose 
width correlates with the strength of the shock flapping mo- 
tions. The fractional widths of the arm, postshock expansion 
zone, and interarm region are typically 10%, 20%, and 70% 
of the arm-to-arm distance, where the gas stays for 15%, 30%, 
and 55% of the arm-to-arm crossing time, respectively. The 
shock flapping motions in the XZ plane make the arm wider 
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than in one-dimensional spiral shocks where the arm takes up 
only 1% of the arm-to-arm distance (Paper II). 

The dense gas produced from TI and shock compression 
tends to sink toward the midplane to form a thin slab, while 
high-altitude regions are dominated by warm rarefied gas. 
The thickness of the dense slab is ~ 10-40 pc, depending 
on the total gas content, heating rate, and presence/absence 
of self-gravity, while the scale height of the rarefied gas is 
Hr 130 pc « crI yjAi^Gpl insensitive to the parameters. 
For model NM.SIO, the thickness ratio of the dense to rar- 
efied components is about 5, which is not much different from 
the results of Dobbs et al. (2008) who reported that the warm 
gas extends vertically up four times more than the cold gas. 
With high density and pressure, the dense slab would trans- 
form to molecular clouds if the appropriate chemical reac- 
tions for molecule formation were included (e.g., Dobbs & 
Bonnell 2007; Dobbs & Price 2008). Thin distributions of 
the cold dense gas are in fact common in numerical simula- 
tions of galactic disks with TI where turbulence is driven by 
magnetorotational instability (Piontek & Ostriker 2007), stel- 
lar feedback via H II regions (Koyama & Ostriker 2009a,b), 
or supernovae explosions (Korpi et al. 1999; de Avillez & 
Berry 2001; Joung & Mac Low 2006; Joung et al. 2009). The 
observed molecular distribution the Milky Way has a scale 
height of ~ 35 pc within the Solar circle, somewhat reduced 
for the most massive clouds (e.g., MaUiotra 1994; Bronfman 
et al. 2000; Stark & Lee 2005). Galactic H II regions are also 
within about 30 pc of the midplane (Lockman 1977). In the 
inner Galaxy and near the Solar circle, the scale height of the 
cold H I layer is about 1 .5 times smaller than the warm H I gas 
(e.g., Falgarone & Lequeux 1973; Crovisier 1978; see also 
Ferriere 2001), although the cold and warm phases appear to 
have a similar scale height in the outer Galaxy (e.g.. Dickey 
et al. 2009). 

4. Random gas motions driven by shock flapping mo- 
tions. — The flapping motions of spiral shocks stir the gas 
and supply random kinetic energy. Allowing for incomplete 
subtraction of streaming motions in the arm region, the in- 
duced density-weighted velocity dispersions are (t^; ^ (jy ~ 
2-3 km s"' for both dense and rarefied components for the 
non-self-gravitating models, with larger values corresponding 
to disks with larger Eo- Compared with the results of Paper 
II where the vertical coordinate was suppressed, these values 
are similar to those for the dense gas in the arms and liirger 
by a factor of ~ 2 - 3 for the interarm rarefied gas. This im- 
pUes that it is the rarefied gas that is more efficiently stirred by 
the shock flapping motions. The self-gravitating models have 
larger velocity dispersions, in the range ax ^ ay ^4-5 kms"' 
for the dense and cr;^ ~ cr^, ~ 4-7 kms~' for the rarefied gas, 
indicating that self-gravity enhances shock flapping and ve- 
locity dispersions, especially for rarefied gas. These in-plane 
velocity dispersions in the current multiphase models are sim- 
ilar to those in the isothermal models considered in Paper I. 

The vertical velocity dispersions of the rarefied gas in NU 
and NM models are cr^ ~ 1 .7 kms~\ insensitive to So- In NU 
models, the vertical motions of the dense gas are excited pref- 
erentially by vertical motions of the rarefied gas. Since the 
mass fraction of the rarefied gas decreases with Eq, the verti- 
cal velocity dispersions of the dense gas in NU models is a de- 
creasing function of Eq, varying roughly as a^ oc S^-^. In NM 
models, the postshock gas is overpressured due to enhanced 
heating and thus expands verticaUy, increasing a^ compared 
to NU models. 



The level of random gas motions in our models are gen- 
erally smaller than the observed velocity dispersions ^ 7 — 
lOkms"' for atomic gas in the solar neighborhood (e.g., 
Heiles & Troland 2003) and for the larger molecular clouds 
in the Milky Way (e.g.. Stark & Brand 1989; Solomon 
et al. 1987; Heyer et al. 2009). Thus, we conclude that 
other sources of the interstellar turbulence (e.g., Mac Low & 
Klessen 2004; Eknegreen & Scalo 2004) must exceed that 
provided by spiral shocks. One of the dominant mechanisms 
is of course supernova explosions (e.g., Korpi et al. 1999; 
de Avillez & Breitschwerdt 2005; Joung & Mac Low 2006; 
Joung et al. 2009). In outer regions of galaxies where star for- 
mation activity is low, the magnetorotational instability (e.g., 
Sellwood & Balbus 1999; Piontek & Ostriker 2005, 2007) 
and/or cosmic infall of gas (e.g., Santillan et al. 2007, 2009; 
Klessen & Hennebelle 2009) may play an important role in 
driving the ISM turbulence. H II region expansion and radia- 
tion pressure are important in injecting energy into the ISM as 
GMCs are dispersed (e.g., Matzner 2002; Murray et al. 2010). 
At large scales, self-gravitating instability with galactic rota- 
tion and shear can drive turbulence at near-sonic levels (e.g., 
Kim & Ostriker 2002; Wada & Norman 2002; Kim & Ostriker 
2007; Agertz et al. 2009). 

5. Effect of self-gravity and properties of self-gravitating 
clouds. — When self-gravity is included, dense gas in SM 
models with Eq > 5 Mq pc~^ suffers from gravitational in- 
stability, eventually forming two large clouds in each model. 
These are separate in the interarm regions, temporarily merge 
in the arm, and then break up into two pieces in the postshock 
expansion zone. These clouds have a radii ~ 45 -60 pc and 
mass ~ (4-9) X 10^ M© each, and are gravitationally bound 
with a virial parameter of a ^ 2. In our present models, we 
have not attempted to include realistic star formation feed- 
back, but instead increase the heating rate at high density to 
prevent collapse. As a consequence, the main support against 
self-gravity comes from thermal pressure. The mean thermal 
sound speed and internal velocity dispersion of the clouds are 
Cci ^ 3.7-5.1 kms"' and a^ ^ 0.9- 1.8 kms"^ respectively. 
For models with Eq = 2 Mq pc"^, self-gravity is insufficient 
to form bound clouds. Nevertheless, self-gravity increases the 
dense gas fraction by a factor of ~ 2 compared to the non-self- 
gravitating counterpart of this model, which in turn decreases 
the vertical velocity dispersion of the dense gas by a similar 
factor 

Formation of self-gravitating clouds in our two- 
dimensional models requires the production of the dense gas 
due to TI, and then additional shock compression. Although 
our present models do not capture the cloud destruction 
process, bound clouds created inside spiral arms may be 
disrupted before they leave the arms if feedback from star 
formation is sufficiently strong (Shetty & Ostriker 2008; 
Wada 2008). Nevertheless, the presence of high-density, 
self-gravitating clouds in the interarm regions opens an inter- 
esting possibility that the spiral shocks - where the diffuse 
gas is strongly compressed - do not necessarily coincide 
with the regions of highest gas density (in gravitationally- 
bound clouds). For example, Patrikeev et al. (2006) found 
strongly-polarized nonthermal radio emission that may trace 
magnetic arms, detected preferentially upstream of the CO 
arms in the inner disk of the Whirlpool galaxy M5 1 (see also 
e.g., Fletcher et al. 2010). We note, however, that the current 
urmiagnetized models are not yet able to provide clues to 
the relation between gaseous and magnetic arms. It will be 
interesting to see how TI, spiral shocks, and reaUstic star 
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APPENDIX 

TWO-PHASE DISKS IN HYDROSTATIC EQUILIBRIUM 

We consider thermally bistable two-phase disks in which a cold, dense layer with surface density T,d with thickness Ho is 
surrounded by a warm, rarefied medium. Since the scale height of the dense layer is very small compared to that of the rarefied 
gas, we approximate the former as razor-thin (Ho ~ 0). We further assume that the mass fraction of the rarefied gas is small, so 
that its self-gravity is unimportant. Let pRiz) denote the density distribution of the rarefied gas. In the presence of the external 
gravity from a stellar disk of uniform density p*, the condition of vertical hydrostatic equilibrium for the warm gas reads 

4^^=-47rGp*z-27rGSz,sign(z), (Al) 

dz 

where Cr is the isothermal sound speed of the rarefied gas, assumed to be independent of z. Integrating equation (Al) over z, one 
obtains 

2hl 



Pr(z) = PR(0)exp 



z^ + —z 



(A2) 



where = c|/(47rGp,) is the Gaussian scale height which the rarefied gas would have in the absence of self-gravity. The surface 
density, Ex, of the rarefied medium is then given by 

/•OO 

S« = 2/ pR(z)dz = ^t^Gnso), (A3) 

where 

Eng = {2n)'l\pR{0) = -^^= (A4) 



■\/2Gp*c% 



is the surface density without gas self-gravity. 



J"(io) = exp(io)erfc(iy^) (A5) 
is the reduction factor, and 

so^^, (A6) 

measures the strength of gravity due to the dense gaseous slab relative to the external vertical gravity (see, e.g., Kim et al. 2002). 
Note that the results of §3 suggest that when two-phase equilibria are established, the interface between dense and rarefied media 
has a constant pressure Ptians, so that we may take p/{(0) = Pr{Q)/c\ = /a-ans/cl since the dense medium has a very small scale 
height (e.g., Piontek & Ostriker 2007), such that 

E. = %^. (A7) 

The scale height of the rarefied medium is given by 

(l+2*o) 



2 _ _ .2 

R ~ TOO 7T1 ~ "j; 



lZ,PR(z)dz 

where equation (A2) is used. 
The condition of mass conservation requires 



TT J"(so) 



(A8) 



so that 



So = Sz) + Ex (A9) 

2c|p, 



So = —J . (AlO) 



For self-gravitating cases, we fix = 0.056 Mq pc ^ and cr = 7 kms \ and solve equations (A7) and (AlO) iteratively to find 
T,R and Ed = ^^o - ^r as functions of Eq. The resulting values of /« = Es/Eq and Hr are plotted in Figure 3 as dashed lines. For 
non-self-gravitating cases, Er = Eng = f'R(0)/(2Gp*c|)'''^ and Hr = hg corresponding to so = 0, plotted as solid Hues in Figure 3. 

Note that since Eo < Eq, sq < 1 if 7I'GEq/(2c|p«) < 1; for our models with = 0.056 M© pc"^ and Eq < 10 M© pc"^, 
7rGEg/(24p*) < 0.2. When sq < 1, J^{so) « 1, so tiiat//s « ftg 128 pc and Er 2.8 Mg pc'^ for p* = 0.056 Mq pc'^ and 

PtransAB = 2100Kcm-3. 
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